(load "line-search-algorithm.ss")

;;;;;;;;;;;;;;;;;;;;; test 1 ;;;;;;;;;;;;;;;;;;
;; test-function: 2d problem
(define (2D-problem example epsilon c1 c2)
  (lambda (method n)
    ;; n: 迭代次数
    (let ([f (car example)]
          [g (cadr example)]
          [v0 (caddr example)])
      (let* ([solution-stream (method f g v0 epsilon c1 c2)])
        (map (lambda (j)
               (vector->list (stream-ref solution-stream j)))
             (range n))
        ;;(plot (frame-set (xlabel "x") (ylabel "y") (grid)) (linespoints ps))
        ))))

;; example1. Rosenbrock function
;; f(v) = 100 (v2-v1^2)^2 + (1-v1)^2; v=(1,1) is the solution.
(define (rosenbrock)
  (define (f v)
    (+ (expt (- 1 (vector-ref v 0)) 2)
       (* 100 (expt (- (vector-ref v 1) (expt (vector-ref v 0) 2)) 2))))
  (define (g v)
    (let ([aux (- (vector-ref v 1) (expt (vector-ref v 0) 2))])
      (vector (+ (* -2 (- 1 (vector-ref v 0)))
                 (* -400 (vector-ref v 0) aux))
              (* 200 aux))))
  (define v0 (vector -1.2 0))
  (list f g v0))
(display "\n---------------- example 1: rosenbrock ---------------\n")

#|
(let ([result-points (map (2D-problem (rosenbrock) 1e-2 1e-2 2e-2)
                          (list SD-method-stream
                                Subspace/with-SD-method-stream
                                (CG/with-SD-method-stream PR+/beta)
                                Quasi-Newton-SR1-method-stream
                                )
                          (list 20
                                20
                                2;;CG seems not work well
                                20
                                ))])
  (apply plot (frame-set (frame-title "Rosenbrock")
                         (xlabel "x") (ylabel "y") (grid)
                         (xrange -1.5 2) (yrange -0.3 1.5)
                         ;;(save "rosenbrock.pdf")
                         )
         (points (list '(1 1)) (title "x*") (pointtype 6) (pointsize 2))
         (map (lambda (ps name) (linespoints ps (title name)))
              result-points
              (list "SD"
                    "Sub"
                    "CG"
                    "Quasi"
                    ))))
;; Sub 最好, 其次 SD, 之后 Quasi, CG失败?
|#
;;;;;;;;;;;;;;;;;;;;;;;;; test 2 ;;;;;;;;;;;;;;;;;;;
;;;;  test-2 function
(define (test-2 example-fun n epsilon c1 c2)
  (let ([example (example-fun n)])
    (map (lambda (name method)
           (time (printf "\n\n**************** ~a: f(x*) = \n ~a\n"
                         name
                         ((car example)
                          (method (car example) (cadr example) (caddr example)
                                  epsilon c1 c2)))))
         (list "SD"
               "Sub"
               ;;"CG"
               "Quasi"
               )
         (list SD-method
               Subspace/with-SD-method
               ;;(CG/with-SD-method PR+/beta)
               Quasi-Newton-SR1-method
               ))))
;; example1: (rosenbrock)
;;(test-2 (lambda (n) (rosenbrock)) 2  1e-2 1e-2  2e-2)

;; example2: (高维二次型) f(v) = 1/2 v * G * v, with G = diag(l1, l2, ..,ln), li>0
(display "\n---------------- example 2: 高维二次型 ---------------\n")
(define (example2 n)
  (define G-diag (vector-map (lambda (ele) (random 100.))
                             (vector-range n)))
  (define x0 (vector-map (lambda (ele) (- (random 200.) 100))
                            (vector-range n)))
  (define G (diag G-diag))
  (define (f v)
    (* 1/2 (vector-dot v (matrix-*-vector G v))))
  (define (g v)
    (matrix-*-vector G v))
  (list f g x0))
(test-2 example2 45    1e-4    1e-2    2e-2)
;;----------      n   epsilon    c1      c2
;; Sub 最好, 其次CG 和 Quasi, 最差 SD

;;; example3: f(v) = Sum_i {e^vi - vi}
(display "\n---------------- example 3: covex-1 function ---------------\n")
(define (example3 n)
  (define (f v)
    (apply + (vector->list (vector-map (lambda (x) (- (exp x) x)) v))))
  (define (g v)
    (vector-map (lambda (x) (- (exp x) 1)) v))
  (define v-init
    (vector-map (lambda (i) (/ i n)) (vector-range 1 (+ n 1))))
  (list f g v-init))
(test-2 example3 1000    1e-5    1e-3    2e-2)
;;                n   epsilon    c1      c2
;; Note: SD 最好, 其次CG,Subspace, Quasi最差.

;;; example4: f(v) = Sum_i {i/10 (e^xi -xi)}
(display "\n --------------- example 4 : covex-2 funtion--------------\n")
(define (example4 n)
  (define (f v)
    (apply + (vector->list (vector-map (lambda (i x) (* i 1/10 (- (exp x) x)))
                                       (vector-range 1 (+ n 1))
                                       v))))
  (define (g v)
    (vector-map (lambda (i x)
                  (* i 1/10 (- (exp x) 1)))
                (vector-range 1 (+ n 1))
                v))
  (define v-init (make-vector n 1.))
  (list f g v-init))
(test-2 example4 1000    1e-3    1e-3    2e-2)
;;----------      n   epsilon    c1      c2
;; Sub,CG 最好, 其次 SD, 最后 Quasi

;;; example5: penalty function, f(v) = a Sum_i (xi - 1)^2 + (Sum_i xi^2 - 1/4)^2, a = 1e-5
(display "\n --------------- example 5 : penalty function --------------\n")
(define (example5 n)
  (define a 1e-6)
  (define (aux-num v)
    (- (vector-dot v v) 1/4))
  (define (aux-vec v)
    (vector-map (lambda (vi) (- vi 1)) v))
  (define (f v)
    (let ([aux-v (aux-vec v)])
      (+ (* a (vector-dot aux-v aux-v))
         (expt (aux-num v) 2))))
  (define (g v)
    (vector-add (scalar-*-vector (* 2 a) (aux-vec v))
                (scalar-*-vector (* 4 (aux-num v)) v)))
  (define v-init (let ([aux (make-vector n 0)])
                   (vector-set! aux 0 (+ 0.5 (random (+ n n 0.0))))
                   aux))
  (list f g v-init))
;;(test-2 example5 20    1e-7    1e-3    2e-2)
;;----------      n   epsilon    c1      c2
;; Sub最好, 其次SD; CG, Quasi失败?

;;; example6: f(v) = Sum_i {n - Sum_j {cos(xj) + i (1-cos(xi)) - sin(xi)}}^2
(display "\n --------------- example 6: 三角函数 --------------\n")
(define (example6 n)
  (define (aux-vec1 v)
    (vector-map (lambda (j vj)
                  (- n (apply + (vector->list (vector-map cos v)))
                     (* j (- (cos vj) 1))
                     (sin vj)))
                (vector-range 1 (+ n 1))
                v))
  (define (aux-vec2 v)
    (vector-map (lambda (j vj)
                  (- (* j (sin vj)) (cos vj)))
                (vector-range 1 (+ n 1))
                v))
  (define (f v)
    (let ([s (aux-vec1 v)])
      (vector-dot s s)))
  (define (g v)
    (vector-add (scalar-*-vector (* 2 (apply + (vector->list (aux-vec1 v))))
                                 (vector-map sin v))
                (vector-map (lambda (a b)
                              (* 2 a b))
                            (aux-vec1 v) (aux-vec2 v))))
  (define v-init (make-vector n (/ 1 n)));;维数越大,初值越小
  (list f g v-init))
(test-2 example6 300    1e-4    1e-3    2e-3)
;;----------      n   epsilon    c1      c2
;;Quasi最好, 其次SD, 之后Sub, CG失败?

(exit)
